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ABSTRACT 

Using data taken as part of the Bluedisk project we study the connection between 
neutral hydrogen (Hi ) in the environment of spiral galaxies and that in the galaxies 
themselves. We measure the total Hi mass present in the environment in a statistical 
way by studying the distribution of noise peaks in the Hi data cubes obtained for 
40 galaxies observed with WSRT. We find that galaxies whose Hi mass fraction is 
high relative to standard scaling relations have an excess Hi mass in the surrounding 
environment as well. Gas in the environment consists of gas clumps which are individ¬ 
ually below the detection limit of our Hi data. These clumps may be hosted by small 
satellite galaxies and/or be the high-density peaks of a more diffuse gas distribution 
in the inter-galactic medium. We interpret this result as an indication for a picture in 
which the Hi-rich central galaxies accrete gas from an extended gas reservoir present 
in their environment. 
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1 INTRODUCTION 

Gas accretion is an important process in the evolution of 
galaxies. Galaxies like the Milky Way would run out of their 
gas in less than 1 Gyr given their current star formation rate 
unless their gas reservoir is continuously replenished (Ken- 
nicutt 1983). The source of this fresh gas is unknown but 
the majority is unlikely to be contained in gas-rich satellite 
galaxies because they can contribute at most 10 percent of 
the required gas mass (Kauffmann et al. 2010; Di Teodoro 
& Fraternali 2014, ). Likewise, high-velocity Hi clouds con¬ 
tribute 3-5 times less gas than required to sustain the Milky 
Way star formation (Richter 2012). A potentially more sub¬ 
stantial source of fresh material is the cooling flow of low- 
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metallicity, ionised gas in a galaxy’s halo (Shull et al. 2009). 
This may be related to the ’’fountain” gas which is driven to 
a height of a few kpc from the disc by star formation, and 
may trigger a gas inflow rate very close to the star forma¬ 
tion rate (Fraternali et al. 2013). Recent theoretical studies 
suggest that the fountain driven gas naturally produces an 
outside-in shrinking galactic disk (Elmegreen et al. 2014), 
which is consistent with what is observed for nearby dwarf 
galaxies (Zhang et al. 2012). However, it does not explain 
the observed inside-out growth of more massive disk galaxies 
(J. Wang et al. 2011), or their chemical evolution (Spitoni 
et al. 2013), which have been both explained in terms of a 
cosmological radial gas accretion (Kauffmann 1996, Spitoni 
et al. 2013). 

In the framework of the standard ACDM cosmology, 
two modes of (radial) gas accretion occur. Gas entering the 
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potential well of a galaxy gets shocked to the virial temper¬ 
ature of the dark matter halo at an early epoch of infall, 
and then gradually cools and falls onto the central galaxy 
(Rees &; Ostriker 1977; Silk 1977; Binney 1977, White &Rees 
1978). Gas that has never been shock heated close to virial 
temperature falls onto the galaxy from outside the virial ra¬ 
dius along filamentary structures (e.g. Keres et al. 2005). 
This latter mode is termed “cold mode” accretion, and is 
believed to be the dominant way of gas accretion for galax¬ 
ies in low mass dark matter halos (Keres et al. 2005, Dekel 
& Birnboim et al. 2006). Numerical studies characterise the 
filaments of cold mode accretion as cool and clumpy clouds 
surrounded by ionised gas (Keres & Hernquist 2009). 

There have been many observational attempts to trace 
the cold mode accreting gas around galaxies. Most of the 
evidence comes from observing ionised gas at intermediate 
to high redshift. QSO absorption lines in Mg II with an 
inflowing velocity feature have been found to be prevalent 
at a distance of tens of kpc around star forming galaxies 
(e.g. Giavalisco 2011, Kacprzak et al. 2012). At least one of 
them is observed to have a filamentary structure (Rauch et 
al. 2011). Furthermore, metal-poor Lyman limit systems are 
observed to have characteristics expected for cold accretion 
(Lehner et al. 2013). 

Tracing the cold-mode accreting gas in the neutral 
phase in galaxies at low redshift has been difficult. Hi clouds 
that form a filamentary structure have been found around 
some nearby galaxies (e.g., Shostak & Skillman 1989; Oost- 
erloo et al. 2007; de Blok et al. 2014, for a review see Sancisi 
et al. 2008) and in the Milky Way (HVCs, Putman et al. 
2002) but are more likely to trace interaction with compan¬ 
ions and star formation feedback. None of them has been 
confidently demonstrated to be tracing cold-mode gas accre¬ 
tion. Although cold-mode accretion is believed to produce 
the highly clumpy galactic disks observed at high redshift 
(Agertz et al. 2009), we do not find an especially clumpy 
morphology for the Hi-rich galaxies when compared to con¬ 
trol galaxies at low redshift (J. Wang et al. 2013, 2014). 
This may be because at low redshift the accreting gas gets 
slowed down and blends better with the hot halo gas around 
the galaxies (Putman et al. 2012). Hence, it may be better 
to search for the cold-mode accreting neutral gas well out¬ 
side the galactic disks, by searching for signal from a large 
volume around the central galaxies. Here we attempt to per¬ 
form exactly this search by analysing the data taken as part 
of the Bluedisk project (J. Wang et al. 2012) in a special 
way. 

Previous analysis of the same Bluedisk dataset already 
shows that the gas richness of a galaxy is linked to the gas 
richness of its environment. In particular, the satellites of ab¬ 
normally Hi-rich spirals are themselves abnormally Hi rich 
(E. Wang et al. 2015). This result may be seen as a new 
manifestation of the “conformity phenomenon”, whereby the 
colours of satellites correlate with the colour of central galax¬ 
ies (Weinmann et al. 2006, Kauffmann et al. 2010, W. Wang 
et al. 2012). The work of E. Wang et al. (2015) was based 
on individual Hi detections in the Bluedisk data cubes and, 
therefore, is limited to galaxies with Mhi > 10 8 M©. As 
previously pointed out by Kauffmann et al. (2010), because 
the galactic conformity usually extends to scales far beyond 
dark matter halos, it implies that the satellite galaxies are 
probably tracing the underlying smooth gas reservoir that 


is available for fuelling the central galaxies. If this hypothe¬ 
sis is true, the conformity behaviour should still hold when 
it is examined at low luminosity levels (i.e. by tracing the 
less bright satellites or clouds in the general environment 
that are potentially far more numerous than the brighter 
satellites). 

Motivated by these results we perform a new analysis 
of the Bluedisk data. Our aim is to trace the faint gas distri¬ 
bution that may surround galaxies below the detection limit 
of the Hi data. Our method consists of adding up the flux of 
(both positive and negative) noise peaks in the Hi data cubes 
out to a radius of 16 arc min (typically 500 kpc), trying to 
detect a statistical Hi signal (i.e. an excess of positive de¬ 
tections) and studying whether its properties correlate with 
the Hi content of the central galaxy. 

The paper is organised as follows. In section 2, we de¬ 
scribe the Bluedisk data used in this analysis. In section 3, 
we introduce our new method to measure a total Hi mass in 
faint systems around primary galaxies through examining 
the noise peaks. In section 4, we discuss the origin of the 
Hi mass measured from an analysis of the noise peaks. We 
find that the excess Hi mass fraction in the primary galax¬ 
ies is correlated with an excess Hi mass in the surrounding 
environment. In section 5, we discuss the implication of our 
results for cold mode gas accretion in galaxies. A ACDM 
cosmology with — 0.3, — 0.7 and h — 0.7 is assumed 

throughout the paper. 

2 THE SAMPLE 

2.1 The Bluedisk dataset 

The Bluedisk project aims at searching for clues about gas 
accretion and galactic disk formation in the local universe. 
Details about the sample, data processing and analysis can 
be found in Paper I and II, and here we only review the 
most relevant information. 

The sample consists of 50 massive galaxies (M* > 
10 lo M©) with a redshift between 0.018-0.03 (correspond¬ 
ing to a luminosity distance of 70-130 Mpc). From the at¬ 
las in Paper I, most of them are blue spiral galaxies. The 
target fields are covered by the Sloan Digital Sky Survey 
(SDSS, Abazajian et al. 2009), the Galaxy Evolution Ex¬ 
plorer (GALEX) imaging survey (Martin et al. 2005), and 
the Wide-Field Infrared Survey Explorer (WISE, Wright et 
al. 2010), allowing a multi-wavelength analysis. Radio 21 cm 
synthesis data were obtained with the Westerbork Synthesis 
Radio Telescope (WSRT), with the main goal of observing 
the Hi emission line. In this paper, we make use of the Hi 
data cubes optimised for high sensitivity at lower resolu¬ 
tion, which are produced with Robust 0 weighting and 30” 
tapering. The data cubes have a typical rms of 0.37 mJy 
beam -1 , a typical beam size of 38^ x 36^ and a velocity res¬ 
olution of 12.3 km/s (after Hanning smoothing), which pro¬ 
vides sufficient sensitivity to detect point sources down to an 
Hi mass of a few 10 8 M©. The primary beam has a FWHM 
(full-width-at-half-maximum) of ~ 0.5°, corresponding to a 
~ 1 Mpc sky region around the target galaxies, hence the 
data is very suitable for studying the relatively large-scale 
Hi environment around the central sources. The details of 
the primary beam correction are discussed in E.Wang et 
al.(2015). 
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Following Paper I, we use the 42 isolated target galaxies 
for this study. Hence, the local environment is different from 
groups or clusters, where the dominant process is gas strip¬ 
ping rather than accretion. We further exclude two galaxies 
(galaxy 14 and 27) with the lowest redshift, so that we have 
a narrow redshift range of 0.023-0.03 (a luminosity distance 
of 100-130 Mpc). We call these 40 target galaxies “primary 
galaxies”. 

These galaxies have a broad distribution in Hi mass 
fraction /hi = logMm/M*, ranging from lower than -2 to 
higher than 1 (Figure [TJ . Compared to other samples fea¬ 
turing Hi-rich galaxies (like the samples from Lemonias et 
al. 2014 and the HIghMass project, High Hi Mass, Huang 
et al. 2014), the Bluedisk sample has the advantage of be¬ 
ing able to compare the most Hi rich galaxies with those of 
normal Hi content. As we will show later in Section T4.21 this 
is crucial for producing our main results. 

As mentioned before, the aim of the paper is to investi¬ 
gate the connection between the Hi in the (central) primary 
galaxies and the Hi in the environment. In addition to /hi, 
we also use a few other parameters to quantify the Hi rich¬ 
ness of primary galaxies. Using the GALEX Arecibo SDSS 
Survey (GASS, Catinella et al. 2010), Catinella et al. (2013, 
C13 for short hereafter) calibrated a photometric esti¬ 
mator of /hi in galaxies with NUV —r <4.5 (which is also 
the colour range of the Bluedisk sample). The estimator is a 
combination of the NUV —r colour and the mass surface den¬ 
sity, reflecting the connection between star formation rate 
and cold gas content. The difference between the observed 
and the estimated /hi, which we denote as Aci3/hi, mea¬ 
sures the excess gas, and can be used to indicate Hi richness. 
Xiao et al. (in prep) further add the optical concentration 
parameter and stellar mass to the estimator, giving more 
control on the internal structure of galaxies. We also use the 
difference from this 4-parameter (4p) estimator ( A^/hi) as 
a measure of Hi richness. 


3 METHOD 

As discussed in the introduction section, the conformity be¬ 
haviour of central and bright satellite galaxies can be inter¬ 
preted in terms of an underlying extended reservoir of cold 
gas. If this hypothesis is correct, we expect to detect confor¬ 
mity also at lower Hi mass levels than probed by E. Wang 
et al. (2015), in a regime dominated by very faint satellite 
galaxies and/or gas clouds in the inter-galactic medium. In 
this section we introduce our technique to detect this faint 
Hi gas, in the environment of the Bluedisk galaxies. 

The method entails an analysis of the noise peaks (both 
positive and negative) in the cleaned data cubes in or¬ 
der to obtain a statistical detection. Because the noise is 
distributed symmetrically around 0 (both the median and 
mean background of a cube are typically within ±10 -5 mJy 
beam -1 , while the rms of a cube is ~0.4 mJy beam -1 ), a 
positive noise peak will statistically be cancelled out by a 
negative noise peak of the same level. If in addition to noise 
there is a collection of very faint Hi objects, then these will 
add a slight positive signal to the noise, that will elevate 
the positive noise peaks or diminish negatiive noise peaks. 
This bias is detectable by a careful analysis of the distribu¬ 
tion of the noise values throughout the cube. We do this by 
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Figure 1. Hi mass fraction of the Bluedisk sample (dark grey 
dots), the HIghMass sample (triangles), and the sample from 
Lemonias et al. (2014) (L14, diamonds). The dashed line shows 
the median relation between /hi and stellar mass for general 
galaxies, taken from Catinella et al. (2010). The dotted line is 
the dashed line 0.6 dex upward, showing the original division be¬ 
tween the Bluedisk high /hi and control galaxies. 


selecting all (positive and negative) noise peaks above 4 a 
in the area of interest, after removing the signal that has 
been reliably detected in the source finding procedure ( see 
section EU). 

Selecting noise peaks above a certain threshold has an 
important advantage compared to integrating over all vox¬ 
els within a large area. Bright HI emission as well as diffuse 
faint HI emission in the data cubes is cleaned by means of 
clean masks (Paper I). The clean masks do not include the 
faint noise peaks used for our statistical analysis here, and 
these are therefore not cleaned. Because the integral over 
the dirty beam over a large area goes to zero it becomes 
very difficult to determine fluxes over extended areas with¬ 
out proper deconvolution. This is inherent to the sidelobe 
structure of the PSF. Any blind measurement of the total 
flux over a large area is therefore fraught with problems. 
However, by pre-selecting noise peaks, we analyse only the 
central part of the PSF and exclude any sidelobes. We are 
therefore unable to detect the Hi mass of a very extended, 
smoothly distributed medium, (this is a consequence of using 
an interferometer with limited short baselines), but we are 
able to detect small scale faint emission using the method 
outlined here. 

In the following sections, we describe our procedure in 
details. We improve the data quality by flattening the spec¬ 
tral baselines (section I3.1D . We use the source-finding appli¬ 
cation SoFiA (Serra et al. 2015) to extract all signal exceed¬ 
ing 4 a. After removing all bright, reliable detections (sec¬ 
tion we are left with “noise” peaks which we analyse 
in order to detect a statistical signal. We perform a series of 
careful tests to examine the robustness of the accumulative 
signal against various data reduction artefacts (section 13.3D . 
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Flattening the spectral baselines 


3.1 

In order to improve the quality of the Hi data cubes we fit 
second-order polynomial functions in moving windows along 
each line of sight. This allows us to obtain global functions 
describing the shape of the spectral baselines without spec¬ 
ifying their functional form. In practice, we convolve each 
spectrum with the Savitzky-Golay filter (S-G, Savitzky & 
Golay 1964), an analytical solution for the local fitting pro¬ 
cess. We adopt a S-G filter with a half width of 70 channels 
(840 km/s), wider than the velocity width of all the galaxies 
in our data cubes. This ensures that the smoothed spectrum 
describes the large-scale shape of the baseline only. In Figure 
2 we show the shape of the filter and demonstrate the con¬ 
volution effect on simple functions, including the presence of 
edge effects. Within our analysis range of ±500 km/s around 
the spectrum centre (channels 34 to 114), the S-G filter con¬ 
serves the shape of the functions with a maximum deviation 
of less than 15%. 

Before convolving the Hi spectra with the S-G filter 
we mask channels that belong to one of the Hi detections 
presented in Paper I (see also Sec. 3.2 below). We refer to 
these channels as the initial flagging set, and replace them 
with values interpolated from a first order polynomial fit¬ 
ting to the un-flagged part of spectrum. We calculate cr m , 
the median absolute deviation of the original spectrum from 
the convolved spectrum, and update the mask by including 
all channels that deviate by more than ±3 cr m from the 
convolved spectrum. A flagged channel is replaced by the 
value interpolated from a first order polynomial fit to the 
un-flagged part of the whole spectrum, if it is in the ini¬ 
tial flagging set, or in a flagged region that is wider than 10 
channels and less than 40 channels from one end of the spec¬ 
trum (potential bright sources that lie near the edge of the 
spectrum and cause incorrect continuum over-subtraction 
during the former data reduction). Otherwise the flagged 
channel (which is not likely to belong to a bright source) 
gets replaced by a first order polynomial interpolation from 
nearby channels. We repeat the above steps of flagging and 
convolution until cr m (with typical value of a few times 0.1 
mJy beam -1 ) varies by less than 10 -5 mJy beam -1 , or after 
10 iterations. Figure [3] shows that this procedure models a 
variety of spectral baseline shapes successfully. 

The final step of this procedure consists of combining all 
spectral baseline models into individual cubes. We smooth 
these cubes with a gaussian kernel with a FWHM of 2 pixels 
in the Ra-Dec direction in order to reduce the fluctuation 
between adjacent pixels (this does not significantly modify 
the PSF (which has a typical FWHM of 3.7 pixels). Finally, 
we subtract the resulting spectral baseline model cube from 
the original cube. After this subtraction, the median rms 
of each data cube drops very slightly from 0.367 to 0.363 
mJy beam -1 (and from 0.372 to 0.362 mJy beam -1 for 
the inner channel range 34 to 114). Hence while flattening 
the baselines, the method does not significantly change the 
global noise properties of the cube, In Figure [U we compare 
two data cubes, using the same position-velocity cut, one as 
obtained before the additional continuum subtraction, and 
one as obtained after the additional continuum subtraction. 
While the noise features look similar, the continuum sub¬ 
traction is obviously improved. 



0 20 40 60 80 100 120 140 

channel 


Figure 2. The S-G filter (top row) and examples of the convolu¬ 
tion effect on simple functions. The area under the S-G filter in 
the top row is normalised. From the second to the bottom row, 
the dark grey curves are linear, second order polynomial and third 
order polynomial functions and a flat line with a spike on one side. 
The magenta curves are the result of the grey curves convolved 
with the S-G filter. The horizontal dashed lines mark the y axis 
position of 0, and the vertical lines mark the channels 34 and 114. 
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Figure 3. Examples of spectra before (black lines) and after ap¬ 
plying the S-G filter (magenta lines). The dashed lines mark the 
zero level. 
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Figure 4. An example of an x-z projection from an original 
Hi cube (top) and that after applying the improved continuum 
subtraction discussed in section [3TT] (bottom). The displayed grey 
scale ranges from -3 to 3 a around the mean of the cube. 


old of reliable sources do not dominate the total signal in 
candidates. 

Using these selection criteria, besides the 40 primary 
galaxies, SoFiA detects in total 59,757 sources from the 40 
cubes, of which, using the scheme as described, 137 are clas¬ 
sified as reliable. 42 of the reliable sources lie within a pro¬ 
jected distance ( D pro j ) of 500 kpc and a systematic velocity 
difference ( AV sys ) of 500 km/s around the primary galaxies. 
95% of them have Mm > 1.62 x 10 8 M©. All of them can be 
associated with an optical galaxy down to a r band magni¬ 
tude of 22 by requiring the centre to be offset by less than 
the size of the Hi beam. The connection between primary 
galaxies and these reliable sources (satellite galaxies) is stud¬ 
ied in detail in E. Wang et al. (2015). One of the key findings 
is that the satellites of Hi-excess galaxies (Aci3/hi > 0; see 
Sec. 2) have higher Hi excess than the satellites of Hi-normal 
galaxies (Aci3/hi < 0). Here we aim to test whether such 
a connection between the Hi content of galaxies and that 
of their environment holds also when studying the fainter, 
candidate Hi sources. 


3.2 Source finding and statistical detection of HI 
in the environment 

We extract sources (both real detections as well as noise 
peaks to be used for our statistical analysis) from the im¬ 
proved Hi data cubes 0 using the source detection software 
SoFiA0. We refer to Serra et al. (2015) for a detailed descrip¬ 
tion of the software. Here, we only outline the most relevant 
steps and describe the key parameter choices. 

SoFiA first convolves the data with various smoothing 
kernels and selects voxels with absolute values above a cer¬ 
tain threshold. We use smoothing kernels with a width of 0 
and 3 pixels in the projected sky direction (roughly the size 
of the synthesised beam), and a width of 0, 3 and 5 pixels 
(corresponding to 0, 36 and 60 km/s) in the velocity direc¬ 
tion (smoothing with a kernel with width of 0 is equivalent to 
no smoothing). Each channel of the cube is weighted accord¬ 
ing to the noise, so that noise variations along the velocity 
axis are removed. This step is very useful to suppress the 
detection of weak RFI (radio frequency interference) resid¬ 
uals. We adopt a detection threshold equal to 4 a. We note 
that the clipping is performed on the absolute values of the 
voxels in the filtered data cubes, meaning that our initial 
source list contains sources with negative total flux. These 
parameter settings for SoFiA are chosen after several tests 
and provide a good balance between suppressing the noise 
and gathering sufficient statistics for the following analysis.. 

SoFiA calculates a reliability index R for each source 
with positive total flux (Serra et al.2012, 2015). We chose 
0.01 Jy beam -1 as a minimum total flux and R>0.99 to clas¬ 
sify a positive detection as a reliable source. We will refer 
to the remaining positive and negative detections (the noise 
peaks) as the “candidates” hereafter. While the division be¬ 
tween reliable sources and candidates is not strict, we will 
see in section o that the candidates close to the thresh- 

1 The data cubes are not corrected for primary beam attenua¬ 
tion. The Hi mass of sources are corrected for the primary beam 
attenuation after they are extracted. 

2 https://github.com/SoFiA-Admin/SoFiA 


We find 6,369 candidates that lie within a D pro j of 100 
to 500 kpc and V sys distance of 500 km/s from the primary 
galaxies, and these will serve as the main sample for analysis 
in what follows. We exclude the inner 100 kpc region because 
some of our Hi rich galaxies have large Hi disks extending 
to ~ 100 kpc. The radial range is always confined within 
the full-width-half-power of the primary beam. Most of the 
candidates have Hi masses (absolute values for the negative 
candidates) within the range of 4.0 x 10 6 M© and 1.4 x 10 8 M© 
(5 and 95 percentiles of the distribution). 

The ensemble of all positive and negative candidate de¬ 
tections results in a net positive Hi signal in the environ¬ 
ment around Bluedisk galaxies. We show this by calculating 
for each cube the accumulative Hi mass (Mm.cand) by sum¬ 
ming up all negative and positive candidates in the volume 
defined above. To further improve the signal-to-noise ratio, 
we calculate the average MHi, C and of the 40 cubes. We cal¬ 
culate the error bars as a combination of the variation over 
the data cubes (varl) and the mean variation within a data 
cube (var2). The variation over the data cubes (varl) is cal¬ 
culated through a standard bootstrapping procedure. Tak¬ 
ing the whole sample of 40 cubes for example, we randomly 
select 40 cubes from the sample (repetition allowed) and 
calculate the average MHi, C and- We repeat this step for 1000 
times and obtain a distribution consisting of 1000 values of 
average MHi,cand- varl is calculated as the variation of this 
distribution. The final error is calculated as \/varl + var2. 
We derive an average MHi,cand of ~ 1.1 x 10 8 M© for the 40 
cubes within the whole analysis volume. 

The question is, however, whether and how significantly 
MHi,cand is affected by remaining artefacts in the data. 
Asymmetrically distributed noise can come from weak RFI 
residuals that are not completely removed during data re¬ 
duction, and can also be produced in the continuum sub¬ 
traction and cleaning processes. Hence in the following sub¬ 
sections we undertake a series of tests to identify and quan¬ 
tify these effects, before moving on to the scientific interpre¬ 
tation in section [4] 
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3.3 Influence of synthesis data reduction on the 
mass in candidates 

We begin this series of tests by verifying that the white 
noise in the data does not accumulate to a significant net 
signal. We test this by applying the same source finding 
and Mnpcand measurements presented in section f372l to the 
Stokes Q cubes. Since no polarised signal in the line data is 
expected, the result is expected to be 0 unless indeed a pos¬ 
itive accumulative noise signal is present. Figure [5] confirms 
that Mnpcand of the Stokes Q cubes is indeed close to 0. 

In addition, we visually inspect the data cubes carefully 
in various projections. The well trained eye is able to spot 
and identify systematics in the data. Because the WSRT 
is a linear E-W array, bright RFI residuals show patterns 
of parallel stripes and bright cleaning residuals show pat¬ 
terns of elliptical rings around bright sources on the Ra-Dec 
projections. Improper continuum subtraction results in neg¬ 
ative regions on the Ra-velocity or Dec-velocity projections 
(see Figure We find no obvious cleaning or continuum 
subtraction residuals. We find some very weak underlying 
RFI patterns in a few channels, but none of our Hi candi¬ 
dates are associated with any of those. This is partly because 
the source finder takes into account noise variations within 
cubes (including those caused by residual RFI; section lT2l) 

We continue our search for possible sources of system¬ 
atic error in our measurement of Mnpcand by investigating 
the effect of incorrect cleaning and continuum subtraction. 
Any systematics related to clean residuals (i.e. the pres¬ 
ence of residual side lobe flux in the HI cubes) should 
correlate with a higher cleaned flux or a larger cleaned re¬ 
gion. This opens the possibility for the following tests. 

We calculate the sum of all positive voxels belonging to 
candidates and the absolute sum of all the negative voxels 
belonging to candidates, and we refer to them as f+ C and and 
f—cand • We calculate the total flux in the reliable sources 
and the total flux in the candidates, and refer to these as 
frei and fcand respectively. We calculate f+ C and , f~cand , 
fcand and frei for each channel map within 500 km/s from 
the systemic velocity and 500 kpc from the position of the 
primary galaxy. These are the velocity range and sky re¬ 
gion within which we perform our analysis. The top row of 
Figure [6] shows that there is no correlation between f— C and 
and frei , suggesting that clean residuals do not significantly 
affect the noise (negative voxels) of the data cube. There 
is a very weak relation between /+ can d and / re z, resulting 
in a systematic linear relation of fcand increasing slightly 
as a function of f re i • We argue that part of the relation 
reflects galaxy-environment connections rather than arte¬ 
facts, as we will show later. If we attribute the relation 
fully to clean artefacts around reliable sources, on aver¬ 
age an Hi mass of 0.11xl0 8 M© related to clean artefacts 
is present in each data cube, which is an order of magnitude 
lower than the average Hi mass contained in the candidates 
(Mm,cand ~ 1.1 x 10 8 M©, section 13721) . In the bottom row of 
Figure [6] we perform a similar analysis as in the top row but 
replace f re i with 7V re z, the number of voxels in the reliable 
sources. We get similar trends, and the linear relation be¬ 
tween fcand and N re i suggests that on average a maximum 
positive bias in Hi mass of 0.07xl0 8 M© is associated with 
the cleaning of reliable sources in a data cube. These tests 
demonstrate that the cleaning process may indeed leave a 


small amount of positive residual flux, but it does not sig¬ 
nificantly affect Mnpcand of the candidates. 

Are there faint systematic residuals from continuum 
subtraction that can not be directly caught by eye? We ex¬ 
tract continuum sources from the continuum maps, with a 
3 a detection threshold. These sources have a 5 percentile 
of 0.94 mJy in the flux distribution. We find that 99% of 
the flux in candidates comes from sight lines without de¬ 
tectable continuum sources. The fraction is higher than 95% 
for both positive and negative candidates, if we count them 
separately. Hence continuum subtraction residuals are not 
likely to significantly affect Mnpcand- 

We also find that the flux in candidates is not corre¬ 
lated with radio flux in continuum sources (with a Pearson 
correlation coefficient of 0.06) or the number of continuum 
sources (with a Pearson correlation coefficient of -0.27). We 
divide the galaxy sample equally into two subsets by the 
total flux in continuum sources, maximum flux in contin¬ 
uum sources and number of continuum sources, and com¬ 
pare their Mnpcand (see Figure 0. We find no differences 
for data cubes having different values for the total or max¬ 
imum radio flux present in the form of continuum sources. 
There is a trend that data cubes with fewer radio continuum 
sources have higher Mnpcand- However, the number of back¬ 
ground continuum sources does not correlate with any of the 
galaxy properties discussed in section 14.21 and we conclude 
that this correlation does not affect our main results. 

Finally, we note that cumulatively there is a weak gradi¬ 
ent of Mnpcand on the low and high V sys sides of the primary 
galaxies, which is unlikely to be physical but does not sig¬ 
nificantly affect the main result of this paper. We refer the 
readers to Appendix A for a detailed investigation of this 
issue. 


3.4 A summary of the method 

To summarise, we used SoFiA to extract candidates with 
typical mass of 10 r M©, which is lower than that of reliably 
detected sources but higher than most of the data reduc¬ 
tion artefacts. We have shown that the total Hi mass in 
candidates, Mnpcand, is not sensitive to artefacts produced 
in the continuum subtraction and cleaning processes of our 
data reduction, and potentially traces a reservoir of low sur¬ 
face density Hi gas in data cubes. In the following section we 
study the nature of these candidates and investigate whether 
Afopcand, he., the gas mass contained in the environment of 
Bluedisk galaxies, correlates with properties of the primary 
galaxies. 

4 RESULTS 

4.1 The nature of the HI mass in candidates 

In this section,we attempt to understand better the nature 
of the signal in the cumulative Mnpcand around the Bluedisk 
primary galaxies. From the top row of Figure [8] we can see 
that the average profile of Mnpcand around the 40 primary 
galaxies rises steadily to ~ 1.1 x 10 8 M© from a projected 
distance of 100 kpc out to 500 kpc, and from a systematic 
velocity distance of 0 out to ~ ±300 km/s. In the bottom- 
left panel, the accumulative Mhi, cand is plotted as a func- 
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Figure 5. The average accumulative MHi, C and from stokes Q cubes as a function of projected distance (left) and systematic velocity 
distance (right) to the primary galaxies. The solid black lines show the average profiles around each primary galaxy, and the grey shades 
show the errors. 
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Figure 6. Relations of fluxes in candidates and the reliable sources in channels. f+ C and is the total flux of positive voxels in candidates, 
and f—cand is the absolute total flux of negative voxels in candidates. f re i and N re j are the total flux and voxel number in reliable 
sources. In each panel, the grey scale map shows the number density of data points, the dots with error bars represent the mean values 
in each x-axis bin, and the dotted line shows a linear fit to the dots, with slope and intercept of the fit denoted (in the first and second 
row, respectively) on top of the panel. 


tion of absolute Hi mass of individual candidates, and we 
can see that MHi, C and is not dominated by the most massive 
candidates. The candidates with Mm > 10 8 M® contributes 
only ~ 1/3 of the net signal. We extract the spectra for 
all the candidates, shift them to have the same central sys¬ 
tematic velocity of 7700 km/s, and stack them. The stacked 
spectrum, as displayed in the bottom-right panel of Fig¬ 
ure 0 is a narrow emission line with a gaussian fit a of 
18.9d=5.7km/s. Removing the resolution effect, the line has 
a width of a ~ 14.3 km/s. If we assume rotational systems 
and apply this velocity width to the baryonic Tully-Fisher 
relation (McGaugh et al. 2000), it corresponds to a baryonic 


mass of 4.4=L2.5x 10 7 M©, similar to the typical Mm of a pos¬ 
itive candidate. The candidates could be low mass satellite 
galaxies, Hi clouds, or Hi clumps in low surface density Hi 
disks. 

If part of MHi,cand probes Hi in galaxies, we expect the 
positive candidates to trace galaxies more closely than the 
negative candidates. As a test, we select galaxies with r band 
flux brighter than 20 mag from the SDSS DR7 photometric 
catalog. For each candidate we search for the optical galaxy 
with the smallest projected distance. We compare the distri¬ 
bution of matching distances for the positive candidates, the 
negative candidates and simulated random positions. In the 
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Figure 7. Averaged accumulative MHi ;C and as a function of distance between 100 and 500 kpc from the primary galaxies. In each panel, 
the primary galaxies are evenly divided into two sub-samples by the median value of the parameter denoted in the left-top corner. The 
red colour is for the sub-sample with higher values of the dividing parameter, and the blue colour is for the comparison sub-sample 
with lower values. The solid lines running through coloured shades show the mean MHi, C and with error bars. The green line shows the 
averaged M HI can( j for the whole sample of 40 cubes, to guide the eyes. The black dotted lines show the position of zero Hi mass. The 
top of each panel shows the KS-test probabilities as a function of radius comparing the distribution of surrounding Hi mass in the two 
sub-samples. A KS-test probability of 0.1 means that there is 90% significance that the two distributions are different. 


left panel of Figure [9l the negative candidates closely follow 
the curve for random positions (with a K-S test probability 
of 0.95, meaning that at a 5% significance the null hypothesis 
of the two distributions being drawn from the same parent 
distribution can be rejected). We find that a slightly higher 
fraction of positive candidates can be matched to optical 
galaxies than the negative candidates at small matching dis¬ 
tance (< 50 arsec, ~1.5 times the FWHM of the Hi PSF), 
although a K-S test probability of 0.28 suggests the differ¬ 
ence is not strong. If we further take only the positive candi¬ 
dates with Mhi > 1.25 x 10 r M© (70 percentile in the mass 
distribution of positive candidates; based on the bottom-left 
panel of Figure O these are the candidates that contribute 
nearly ah the net MHi, C and signal) into account, we find a sig¬ 
nificant difference in the distribution of matching distances 
from the negative candidates (K-S test probability ~ 0.04) 
and from the simulated random positions (K-S test proba¬ 
bility ~0.01). We further check if this small excess of optical 
matches for the bright positive candidates might be due to 
residual from the continuum subtraction around continuum 
sources associated with the optical sources. In the middle 
panel of Figure 0 we match the projected position of neg¬ 
ative, positive and bright positive candidates to continuum 
sources. We find the three types of candidates to be indistin¬ 
guishable in the distribution of matching distances when the 
distance is smaller than 50 arcsec. The K-S test probability is 
0.98 for the comparison between negative and positive can¬ 
didates. The K-S test probability is 0.14 for the comparison 
between negative and bright positive candidates, suggesting 
a weak difference, but the curve shows that the bright posi¬ 
tive candidates are less likely to be found along sight lines of 
continuum sources than the negative candidates. Hence we 
can conclude that at least the brightest ~ 1/3 of the positive 
candidate sample appears to be connected to galaxies, but 
not necessarily radio bright galaxies. 

We also expect the positive candidates to have a broader 
line width than the negative candidates if they are more 
likely to be associated with galaxies. This expectation is con¬ 
firmed in the right panel of Figure [9] The K-S test probabil¬ 


ity is 7xl0 -6 for comparison between positive and negative 
candidates, and is 0.04 for comparison between the positive 
and negative candidates with absolute Mhi > 1.25 x 1 O 7 M 0 . 

In summary, by measuring MHi,cand, we add up the sig¬ 
nal from unresolved, low Hi mass objects distributed in a 
large volume (~500 kpc) around the primary galaxies and 
other reliably detected objects. It is very likely that at least 
a fraction of these candidates are small galaxies. The re¬ 
mainder might be galaxies too (e.g., objects below the 20 
mag cut in r band used here), or small Hi clouds in the 
circumgalactic medium. 

4.2 The HI mass in candidates and the properties 
of primary galaxies 

In this section, we compare MHi, C and between galaxies with 
different properties, especially between galaxies that are rich 
in Hi and normal galaxies. We assume that the Hi-rich (ei¬ 
ther with high /hi or Hi excess) galaxies have been accreting 
gas, and hope to find clues about gas accretion by studying 
Mni,cand around the galaxies. 

In Figure flOl we divide the sample of primary galaxies 
evenly into two sub-samples, focusing on different proper¬ 
ties, and investigate their difference in averaged accumu¬ 
lative MHi,cand as a function of Dproj between 100 and 500 
kpc. The dividing galactic properties include the stellar mass 
(M*), the stellar mass surface density (//*), optical concen¬ 
tration (rgo/rso), NUV—r colour, the Hi mass (Mhi), Hi 
mass fraction (/hi), and Hi excess (Aci 3 /hi and A^/hi) 

EL 

The most prominent trend is the sub-sample with lower 
M* having higher MHi, C and than the comparison sub-sample, 
with a higher than 90% significance that the two distribu¬ 
tions are different at a D pro j of 500 kpc. The sub-samples 

1 We point out the pairs of subsamples investigated here all have 
similar ranges of redshift (e.g. the difference in median luminosity 
distance is always less than 6%) 
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Figure 8. Top row and bottom-left panel: the average accumulative MHi, C and as a function of projected distance and systematic velocity 
distance to the primary galaxies, and as a function of the individual Hi mass (absolute value) of candidates. The solid black lines show 
the average and the grey shades show the errors. Bottom-right panel: stacked spectrum of candidates (black line) and its gaussian fit 
(red line). 
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Figure 9. The onnection of negative (orange), positive (green) and bright positive (purple) candidates with optical and radio continuum 
sources. Left panel: accumulative fraction of projected distance between the candidates and their closest optical galaxies. The dashed line 
is built from matching simulated random positions with the optical galaxies. Middle panel: accumulative fraction of projected distance 
between the candidates and radio continuum sources. Right panel: accumulative fraction of velocity width (length of the detection mask 
in the velocity direction). The colours used for the different sub-samples are denoted in the left-top corner. 
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with higher /hi and A^/hi also have considerable higher 
MHi,cand than the comparison sub-samples; especially the 
latter has a close to 90% significance for the difference in dis¬ 
tribution from the comparison sample. There are also weak 
trends that the sub-sample with lower /x* and bluer NUV—r 
colour have higher MHi, C and (considering the difference in 
average values and the close-to 0.2 K-S test probabilities). 
There is no significant correlation between Mnpcand and the 
optical concentration, the Hi mass or Aci 3 /hi of primary 
galaxies. 

We further control stellar mass to reveal a more intrin¬ 
sic correlation between Mni, C and and Hi richness of primary 
galaxies. For each of the sample dividing parameters, we 
derive the linear relation with stellar mass, and calculate 
vertical distance to the relation. We use this distance to di¬ 
vide the primary galaxies evenly into two samples, and study 
their difference in averaged accumulative Mnpcand- In Fig¬ 
ure [Tl] we only consider the dividing parameters which have 
already shown an indication for a correlation with Mni,cand 
in Figure [10] We can see the trend of galaxies with higher 
/hi having on average higher Mni, C and becomes more promi¬ 
nent as suggested by the K-S test probabilities. The trend 
as a function of A^/hi remains almost unchanged, suggest¬ 
ing the advantage of A^/hi as an unbiased measure of Hi 
richness. Bluer primary galaxies also have on average higher 
Tfnpcand- Thn,cand is no longer dependent on /x* of primary 
galaxies. These trends are consistent with each other, for at 
the same stellar mass, Hi-rich galaxies are also more star 
forming galaxies. Finally, we note that, when the sample is 
divided by /hi, A^/hi and NUV—r of the primary galaxies, 
the difference in the Mnpcand distribution of the sub-samples 
is tentative, with a K-S test probability of slightly above 0.1 
at a Dproj of 300-500 kpc. 


5 DISCUSSION: EVIDENCE OF GAS 

CONFORMITY 

We have found evidence for the presence of an excess of 
faint gas clouds in a volume out to 500 kpc in projected 
distance and 500 km/s in velocity around Hi excess galax¬ 
ies, i.e., galaxies that have more Hi mass than other galaxies 
with similar optical properties. A parallel study of individual 
satellite galaxies around the Bluedisk primary galaxies re¬ 
vealed that satellites around Hi excess galaxies have higher 
Hi excess than satellites around galaxies with normal Hi 
contents (E. Wang et al. 2015). Although these two studies 
are based on different signal-to-noise levels of the Hi data, 
reflecting different reservoirs of cold gas, they consistently 
tell the same story: the Hi reservoir in the environment of 
(massive) galaxies appears to reflect the Hi content of those 
galaxies. As discussed in K10, this environmental Hi may 
trace the dense gas tunnelled into halos through filaments, 
the so-called “cold mode” accretion (Keres et al. 2005). 

Recent numerical studies have shown that cold mode ac¬ 
cretion is not affected by galactic winds (Powell et al. 2011), 
hence it can easily leave its footprint on the low-mass satel¬ 
lite galaxies. Beside the gas locked into satellites, our Hi 
candidates may also include the Hi clouds that condense 
from the ionised cold mode gas (Keres et al. 2009a, 2009b). 
It is interesting to note that the lower limit in the Hi mass 
of candidates detected in our data (~ 10 6 M 0 , imposed by 


our 4-cr detection threshold) is roughly the lower mass limit 
required by an Hi cloud to survive in a hot gas halo around 
galaxies (Murray &Lin 2004). Due to the lack of baseline 
spacings short-ward of 36 m, we are unable to detect the 
very extended, smoothly distributed gas. However, the de¬ 
tected candidates that trace the peaks in the IGM or small 
companion galaxies, indicate that an underlying reservoir of 
gas from cold mode accretion might be present. 

In the picture of cold mode accretion, the outer region 
of the galactic disks is the favoured location of the accre¬ 
tion taking place (as in the classical model of hot accretion) 
(Pichon et al. 2011). The accretion rate is bursty because 
the inflow filament is clumpy (Brooks et al. 2009). These 
two factors may produce the excess Hi in Hi-excess galaxies. 
The Bluedisk primary galaxies have stellar masses around 
that of the MW, for which around 50% of the gas accre¬ 
tion occurs in the cold mode (van de Voort & Schaye 2012). 
The most prominent trend found in this paper is that the 
primary galaxies with lower stellar masses have higher Hi 
mass in candidates around them, which may suggest a tran¬ 
sition between cold- and hot-mode accretion, as predicted 
by the simulations (Dekel & Birnboim et al. 2006). Recently 
Kauffmann (2015) found a tendency for the satellites to 
align along the major axis of late-type Hi-rich central galax¬ 
ies with M* < 10 10-5 M©; they found no such phenomenon 
for more massive central galaxies. Their result also suggests 
such a transition. The median M* of our primary galaxies 
(1 q10.56^q ^ j g consistent with their suggested transition 
point. 

An alternative picture is that the accumulative signal 
from Hi candidates traces a rather smooth accretion of dark 
galaxies. Genel et al. (2010) studied the growth of dark mat¬ 
ter halos and found that at least 40% of the baryons in 
the halos are accreted from very low mass halos, contain¬ 
ing smooth cold T ~ 10 4 K gas and no stars. Although no 
dark galaxies have been found so far from nearby surveys 
(Zwaan& Briggs 2000, Zwaan 2001, Koribalski et al. 2004, 
Haynes 2007), this can be explained if dark galaxies are 
preferentially found around gas accreting massive galaxies: 
galaxies with very high /hi and as massive as the Bluedisk 
sample are rare in the nearby universe. This picture may still 
be consistent with the cold accretion scenario, since recent 
simulations show that halos in filaments may have more sub 
haloes than those in other environment (Guo et al. 2014). 

To this date no gas accreting filaments or clouds have di¬ 
rectly been observed to exist within a few tens of kpc around 
galaxies in the local universe, even when the observational 
depth reaches 10 19 atom cm -2 , as in the case of the WSRT 
Hydrogen Accretion in LOcal GAlaxies Survey (HALOGAS, 
Heald et al. 2011). However, numerical simulations demon¬ 
strate that the accreting gas may just get slowed down and 
ionised through interaction with the hot gas halo when it 
gets close to the galactic disks (Putman et al. 2012). In our 
study here, the candidates are beyond 100 kpc from the 
galactic centre, far above the disk-halo interface. 

Finally, we address the question if the merging of satel¬ 
lites contributes significantly to the cold gas in primary 
galaxies. We select the 20 primary galaxies with higher than 
median /hi, and only 8 of them have reliably detected satel¬ 
lites within a D pro j of 500 kpc and AV sys of 500 km/s. We 
calculate the timescale for a pair of galaxies to merge follow- 
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Figure 10. Averaged accumulative Mni^and as a function of distance between 100 and 500 kpc from the primary galaxies. In each 
panel, the primary galaxies are evenly divided into two sub-samples by the median value of the parameter denoted in the left-top corner. 
The red colour is for the sub-sample with higher values of the dividing parameter, and the blue colour is for the comparison sub-sample 
with lower values. The solid lines running through coloured shades show the mean MHi, C and with error bars. The green solid lines show 
the averaged accumulative Mni^and °f the whole sample for guiding the eye. The black dotted lines show the position of zero Hi mass. 
The top of each panel shows the KS-test probabilities as a function of radius comparing the distribution of surrounding Hi mass in the 
two sub-samples. The dotted lines show the positions of 0, 0.1 and 0.2. A KS-test probability of 0.1 means that there is 90% significance 
that the two distributions are different. 


ing the empirical formula from Kitzbichler & White (2008), 


T = 1.6 


M* 


25h~ 1 kpc K 3 x 10 10 /i —1 A/© 


)-°- 3 Gyr (1) 


where M* is the total stellar mass of the pair. We replace M* 
in the formula with Mhi plus M* to minimize the timescale. 
We calculate the gas merging rate, and compare it with the 
SFR for these 8 primary galaxies with satellites. The highest 
ratio is 0.12, hence the merging rate is too low to support 
the SFR, and can not be a major source for replenishing the 
gas storage in the Bluedisk Hi-rich galaxies. It is consistent 
with the findings of many previous studies like K10 and Di 
Teodoro &Fraternali (2014). 

The gas-rich satellites studied by K10 and E. Wang et 
al. (2015) might be the tip of the iceberg of the missing gas 
accretion. Here with the detection of an excess in the low- 
mass low-column density Hi environment, we might have 
detected a substantial part of the iceberg below the sea- 
level. 


6 SUMMARY AND FUTURE PROSPECTS 

We have developed a new technique to extract information 
about the large-scale distribution of Hi mass from 21 cm 


synthesis data cubes through a technique similar to stack¬ 
ing, which works in an Hi mass regime which is lower than 
achievable for reliable direct detections, but still high enough 
not to be dominated by data reduction artefacts.. We used 
our technique to investigate the mass contained in candi¬ 
date Hi sources, possibly associated with individual satellite 
galaxies, around the Bluedisk primary galaxies. We found 
a connection between Hi excess (A^/hi, or high values of 
/hi at a fixed stellar mass) in primary galaxies and excess 
of Hi in low mass candidates in the surroundings. These Hi 
masses may be the tip of the iceberg of an underlying ex¬ 
tended reservoir of gas that fuels the primary galaxies. It is a 
direct detection of the galactic gas conformity phenomenon 
down to very low Hi masses (~ 10 7 M©). The result is con¬ 
sistent with the cold mode accretion in cosmological simula¬ 
tions. In such a picture, primary galaxies and satellites are 
both fuelled by the cosmic web and exhibit the conformity 
phenomenon of gas richness. The unusually blue outer discs 
in Bluedisk galaxies with high Hi mass fractions, indicat¬ 
ing inside-out formation, can also be explained with such 
cosmological gas accretion (Kauffmann 1996). 

At this moment, we are unable to establish general re¬ 
lations between the mass, structure and morphology of the 
stellar disks and the surrounding Hi masses because we 
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Figure 11. The top row demonstrate how the sample is divided with stellar mass controlled (see text). In each panel, the dashed line is 
the linear fit to the whole sample, and the red and blue dots are with vertical distance above and below the median to the dashed line. 
The bottom row is similar as Figure fTOl but the sub-samples are divided with stellar mass controlled as shown in the panel on top. 


are limited by selection effects in our sample of primary 
galaxies. Our sample misses the most massive galaxies (with 
M* > 10 11 Mq) in the universe. In more massive halos, cold 
gas around primary galaxies may behave in a different way 
from what we find here, for those galaxies are more strongly 
affected by the halo (Catinella et al. 2014). Limited by pri¬ 
mary beam effects, we are also unable to investigate how far 
away from the primary galaxies the conformity phenomenon 
extends. Limited by the small sample size, we are unable to 
statistically measure the total Hi mass in all the satellite 
galaxies. Because the accumulative signal (MHi, C and) we ob¬ 
tain is only ~3 a above 0, we are unable to calculate spatial 
gradients, to quantify spatial distribution, or to estimate 
inflow rate. It also remains unclear how the accreting gas 
arrives at the Hi-excess galaxies, and how much of the Hi 
mass in candidates is not locked into optical galaxies (i.e. 
in the form of clumps in the more diffuse regions in and 
between cosmic filaments, or in dark galaxies). 

We note that the emphasis of this paper is the explo¬ 
ration of a new technique. We demonstrate the promising 
potential to investigate Hi properties in the radio synthe¬ 
sis cubes below 5-6 cr, which is generally considered as a 
threshold for detection reliability. We have tried to control 
the systematics as best as possible and reach the tantalis¬ 
ing result that we may have detected evidence for the pres¬ 
ence of excess Hi in the intergalactic medium around Hi 
rich galaxies. We look forward to applying the technique 
and analysis presented in this paper to a much larger, more 
uniform and complete dataset, especially from the ASKAP 
Hi All-Sky Survey, known as WALLABY, and the proposed 
WSRT Northern Sky Hi Survey (WNSHS), both described 
in Koribalski (2012). 
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APPENDIX A: M H i, C and ON THE LOW AND 
HIGH Vsvs SIDES OF THE PRIMARY 
GALAXIES 

In the main part of this paper, we report the detection of 
Hi mass (MHi, C and) by adding the candidates outside the 
galaxies detected in Hi cubes (Sect. 3.3). We find MHi, C and 
to be higher around galaxies with high Hi excess. We have 
performed tests to show that MHi, C and is unlikely to be 
from CLEANing or continuum subtraction residual artefacts 
(Sect. 3.4). However, as shown in the top-left panel of Fig¬ 
ure EH there is a small gradient in Mni, C and between the 
blue and red shifted sides of the Hi cube around the primary 
galaxies. Consequently, part of the difference in MHi,cand 
found between sub-samples defined by properties of primary 
galaxies may be related to this gradient. Indeed, the gradient 
is slightly stronger for the cubes with low A^/hi primary 
galaxies than the cubes with high A 4 P fm primary galaxies 
(the second and third panels in the top row of Figure (All) : 
and the difference of MHi,cand around primary galaxies with 
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Figure A2. Relation between Mni^cand gradients and Modified 
Julian Date (MJD). The vertical dotted line marks the division 
of the year 2011 and 2012. The horizontal dotted line marks the 
position of 0. 


low and high A 4 P /hi is less significant for the low V sys side 
of the cubes than for the high V sys side of the cubes. 

From inspection of the position of candidates in Hi 
cubes, we find that this gradient is not clearly associated 
with any failures of continuum subtraction, CLEANing or 
RFI pattern removal. We do not find such gradients in the 
test with stokes Q cubes (Sect 3.4). However, we find this 
gradient is weakly correlated with the date of observation 
(Figure IA2D with a Pearson correlation coefficient of 0.3, 
mostly driven by the data points for the year 2011, where 
the distribution of gradients is strongly biased toward nega¬ 
tive values. If we exclude the 12 data cubes observed in 2011, 
the gradients of MHi,cand is much weaker for the whole sam¬ 
ple, and for the primary galaxies with low and high A 4 P fm 
as well (the first 3 columns in the bottom row of Figure [All) . 
Meanwhile, the difference in MHi, C and between the two sub¬ 
samples with low and high A^/hi becomes more significant 
for the low V sys side of the cubes, and is similar for the high 
V sys side of the cubes (the fourth and fifth columns in the 
bottom row of Figure [All) . 

Hence we argue that while the cause for the gradient of 
u,cand 1s not found, the effect is not likely to significantly 
affect our main results, Because it is unclear whether the 
gradient effect is an artefact or a statistic coincidence, we 
still use the full sample of 40 data cubes for analysis in the 
paper. Nevertheless, we warn the readers of this uncertainty 
in our experiment and intend to investigate it with larger 
samples in the future. 
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Figure Al. Similar to Figure fTOl The candidates are grouped according to their systematic velocity with respect to the primary galaxies, 
and the A 4 p /hi of the primary galaxies in the cubes. The colors of the shades are for the sub-samples denoted in each panel. The top 
row is for the whole sample of 40 cubes and the bottom row is for the 28 cubes observed in the year 2012. 
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